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Abstract 

The  availability  of  fast  algorithms  for  computing  the  orthogonal  transforms 
such  as  the  fast  Fourier  transform  (FFT),  Walsh-Hadaraard  transform  (WHT) , 
discrete  Chebyshev  transform  (DCT),  and  related  transforms  has  made  such  trans¬ 
forms  increasingly  popular  in  data  processing.  This  paper  is  concerned  with 
the  comparative  evaluation  of  these  transforms  as  applied  to  the  ACDA  teleseismic 
data.  Computer  results  and  program  listings  are  provided.  The  paper  is  concluded 
with  the  recommendation  that  better  orthogonal  transforms  are  needed  for  bot 
signal  processing  and  seismic  discrimination. 
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Comparison  of  Orthogonal  Trans foras  for  Teleseismic  Data 

C.  H.  Chen 

I.  Introduction 

In  recent  years  there  has  been  an  increasing  interest  with  reBpect  to  using  a 
class  of  orthogonal  transforms  in  the  areas  of  pattern  recognition  and  digital 
signal  processing.  In  pattern  recognition,  orthogonal  transforms  enable  a  noninvertible 
transformation  from  the  pattern  space  to  a  reduced  dimensionality  feature  space.  This 
allows  a  classification  scheme  to  be  implemented  with  substantially  less  features, 
with  only  a  small  increase  in  classification  error.  In  digital  processing,  orthogonal 
transforms  implemented  by  fast  algorithms  are  basic  operations  for  digital  filtering, 
spectral  estination,  etc.  In  discrete  (Wiener)  filtering  for  signal  parameter 
estimation,  orthogonal  transforms  "compress"  the  useful  data  to  a  substantially  small 
number  of  elements  and  thus  simplifies  the  filter  structure  and  reduces  the  computation 
load.  Orthogonal  transforms  are  also  useful  in  designing  multiplexing  communication 
systems . 

The  Walsh-Hadamard  transform,  discrete  Fourier  transform,  the  Haar  transform,  the 
slant  transform,  and  discrete  Chebyshev  transform  (also  called  discrete  cosine  trans¬ 
form)  have  been  considered  for  various  applications,  since  these  are  orthogonal  trans¬ 
forms  that  can  be  computed  using  fast  algorithms.  The  Karhunen-Loeve  transform  (KLT) 
is  an  orthogonal  transform  which  is  optimal  in  the  minimum  mean  squared  error  sense; 
but  is  computationally  difficult  .  In  this  paper  we  are  concerned  with  the  comparative 
evaluation  of  these  orthogonal  transforms  for  teleseismic  data  study  including  signal 
processing  and  seismic  discrimination  (pattern  recognition). 

II.  The  Orthogonal  Transforms 

Of  all  the  orthogonal  transforms  presently  available,  the  discrete  Fourier  trans¬ 
form  is  probably  the  most  conveniently  available  and  provides  the  most  familiar  frequency 
domain  informations  on  the  seismic  data.  The  fast  Fourier  transform  also  serves  as 
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the  groundwork  for  all  transforms  with  sinusoidal  basis  functions  such  as  the 
discrete  Chebyshev  transform  and  nonrecursive  digital  filtering.  The  phase 
spectrum  which  is  often  ignored  contains  much  useful  information.  To  utilize  both 
amplitude  and  phase  spectra  in  the  Fourier  domain  is  not  a  simple  task.  The 
homomorphic  deconvolution  or  the  complex  cepstrum  accomplishes  this  objective  by 
taking  the  logarithm  of  the  complex  Fourier  spectrum.  The  problem  with  the  cepstral 
approach  is  that  it  is  difficult  to  obtain  accurate  long  pass  and  short  pass  system 
outputs  especially  with  the  noisy  data.  The  teleseismic  data  is  highly  sinusoidal 
in  nature  and  thus  the  Fourier  domain  transforms  appear  to  be  most  appropriate. 

Once  digitized  for  computer  compatibility,  the  amount  of  data  needed  to  analyze  the 
seismic  waveforms  can  be  reduced  by  the  orthogonal  transforms.  A  fast  Fourier 
transform  reduces  the  data  needed  to  ai^Lyze  the  waveform  by  a  factor  of  10  (it  would 
be  a  factor  of  5  if  both  amplitude  and  phase  spectra  are  used)  with  a  very  small 
percentage  of  the  reconstructed  root-mean-squared  error.  The  fast  Fourier  trans¬ 
form  accomplishes  the  data  compression  by  changing  the  usual  time  domain  represen¬ 
tation  of  the  seismic  signal  to  one  of  frequency  representation.  A  representation 
of  the  seismic  wave  requires  less  data  points  in  the  frequency  domain  than  in  the 
original  time  domain.  A  Walsh-Hadamard  transform  may  also  be  used  in  a  similar 
manner  to  yield  a  data  compression  factor  of  four.  (The  factor  is  reduced  if  the 
plus  and  minus  signs  are  included.) 

Haar  transform  and  slant  transform  are  more  suitable  for  image  processing  than 
for  seismic  study.  The  discrete  Chebyshev  transform  has  the  important  advantage  in 
that  it  is  closest  to  KLT  in  mean  square  error  while  computable  by  fast  algorithms. 
It  does  require  much  more  computation  time  than  Fourier  or  Walsh  transforms .  The 
Walsh  transform  is  computationally  the  best.  The  data  compression  factor  for  DCT 
is  nearly  the  same  as  for  FFT.  The  DCT  should  be  considered  as  a  member  of  Fourier 
domain  transforms  with  good  promise  for  wide  range  of  applications. 
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For  seismic  data  study,  the  minimum  mean  squared  error  criterion,  however,  is 
of  questionable  value.  Such  criterion  depends  only  up  to  the  second  order  statistics 
(the  mean  and  covariance  functions)  which  may  be  not  enough  to  characterize  the 
seismic  waves.  Furthermore,  reconstruction  of  the  seismic  wave  is  not  an  important 
objective.  To  estimate  some  parameters  such  as  to  determine  the  P-wave  amplitude, 
and  to  discriminate  explosion  from  earthquake,  are  among  the  important  objectives 
of  the  seismic  study.  Both  DCT  and  WHT  spread  the  signal  energy  to  higher  spectral 
components  than  the  FFT  and  thus  may  not  be  as  effective  in  dimensionality  reduction 
for  pattern  recognition  as  compared  with  FFT.  Note  that  DCT  is  better  than  WHT  in 
this  sense.  Unfortunately,  the  reduced  feature  space  provided  by  all  orthogonal 
transforms  is  still  too  large  in  dimension  and  the  subsequent  classifier  is  too 
complex  in  structure  or  computation.  For  two-dimensional  display,  the  transformed 
data  has  to A considerably  compressed  further  to  make  the  display  possible. 

Therefore  we  conclude  that  the  presently  available  orthogonal  transforms  have 
to  be  considerably  improved  for  seismic  study  in  processing,  display,  and  dis¬ 
crimination,  At  the  present  time,  the  Fourier  domian  transforms  are  still  most 
useful  for  seismic  study.  New  transforms  must  be  computable  in  fast  algorithms 
while  drastically  compress  the  data  so  that  the  dimension  of  the  reduced  feature 
space  falls  below  10,  which  is  a  factor  of  100  in  reduction  of  1021+  data  points 
typically  used  in  digital  seismic  data  study. 

III.  Computer  Results 

The  seismic  data  considered  is  the  ACDA  Seismic  Signature  Data  Base  provided 
by  Lt.  Colonel  Russell  B.  Ives.  A  complete  spectral  plot  of  all  269  records  was 
submitted  to  him  with  report,  "Fourier  Spectral  Plot  of  ACDA  Seismic  Signature 
Data  Base",  dated  September  12,  1971*.  Some  conclusions  of  the  spectral  analysis 


are: 


1)  With  a  few  exceptions,  all  explosions  have  a  pronounced  spe.ctral  peak 
at  0.5  -  0.7  Hz. 

2)  The  earthquakes  have  many  discrete  frequency  components,  i.e.  they  are 
rich  in  subharmonics. 

3)  The  explosion  spectra  spread  to  h.^ier  frequencies  than  the  • arthquake 

spectra.  ! 

1*)  Most  records  have  strong  DC  components  which  could  be  removed  to  provide 

better  discrimination  results. 

An  important  reason  for  providing  spectral  plot  of  all  records  is  to  enable 
us  to  get  an  accurate  interpretation  of  the  results.  In  the  attachments  of  this 

report  we  have  included  the  following: 

1.  Computer  program  to  plot  the  original  seismic  waveform  with  1024  data 
points  (samples)  at  10  samples  per  second,  and  to  plot  the  amplitude  of  FFT  up  to 
5  Hz. 

2.  Application  note  on  DCT, 

3.  Computer  program  to  plot  WHT  and  DCT. 

4.  A  complete  plot  of  DCT  and  WHT  for  each  record.  On  each  page  there  are 
two  graphs,  the  upper  one  is  for  DCT  up  to  512  data  points  and  the  lower  one  is 
for  WHT  up  to  1024  data  points. 

The  following  are  conclusions  which  c an  be  drawn  from  the  plot  of  DCT  and  WHT: 

1)  DCT  appears  as  a  low-passed  filtered  WHT.  Theoretical  explanation  of  this 
phenomenon  is  not  available.  There  might  be  some  theoretical  relationship  between 
the  two  which  is  not  yet  discovered. 

2)  The  DCT  for  explosion  is  a  more  smooth  function  of  the  number  of  data  points 
than  the  earthquake. 

3)  The  WHT  is  nearly  periodic  (except  in  amplitude)  with  a  period  of 
approximately  256  data  points , 
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DIMENSION  C<257) 

COMPLEX  E(2000> 

_ 

BYTE  IA<  70) ,  IX<9),  IQ<10),  TEST 
EQUIVALENCE  <E,  X) 

DATATiST,  IX ,  10/  '  X  E ',  -X",  'P' ,  'L' »  •  O' S' ,  'Q'r'H't 

♦  '£•  ,  A*  ,  R',  'T  ,  'H'  fk't'K't'E.'/ 

CALL  XYINIT<B,  2000) 

CALL  SIZE(C)  1024 )  _ _ _ 

CALL  CHSIZEIO.  4,  90.  0) 

FORMAT (13) 

IC-0 

READ <6.  6)  IF 
DO  9  1-1,  2000 
X  <2,  l)-0.  O 
X<1,  1  >-6.0 
IC-IC+1 

READ<1,  END-10)  I A 
READ <  1 )  NPT,  (X(hl)i  I-1,NPT> 

IF <  IF.  GE.  IC)  00  TO  8 
REWIND  1 

Gib  TO  7  "  '  . . . .  . .  .  '  ~  ' 

IF<  IF.  NE.  IC)  00  TO  11 
CALL  IXYPT  <230/  0,  0) 

CALL  XYCHAR<  I  A/  20) 

CALL  IXYPT <486,0,  0) 

IF<JA<7_0>.  EG.  TEST)  GO  TO  2  __ 

CALL  XYCHARUG,  10) 

GiO  TO  3 

CALL  XYCHARC IX,  9) 

CALL  IXYP1  <500,  0,  0) 

CALL  IXYPT<21?0,  0,  1) 

XMAX-X< 1,1) 

XMIN-XMAX 
DO  4  1*1,  NPT 

IF < XM1N.  GT.  X<  1,  I ) )  XMIN*X <  1 ,  I ) 

IF<  XMAX.  LT.  X<  1,  I ) )  XMAX*X  <  1,  I ) 

XS-1690.  0/ <XMAX-XMIN) 

N— 2 190.  0+XMIN#XS 
CALL  IXYPT <N,  0,  0) 

CALL  IXYPT<N,  4000,  1) 

IXX-0 

I  Y*2190.  0+<  XMIN-X  <1,  1 )  )*XS 
CALL  IXYPT <  IY,  IXX,  0) 

DO  5  1-2,  NPT 

TxX-I  +  I-2  . . "  '  '  ■  ■  . 

IY-2190.  0+<  XMIN-X  <  1,  I)  )#XS 
CALL  IXYPT<IY,  IXX,  1) 

CALL  IXYPT <2390,  0,0) 

CALL  IXYPT <4090,  0,  1 ) 

CALL  IXYPT < 4090,  4095,  1 ) 

‘  NP-NPT-1024  "  “  . . ~~~ ' 

IF<NP.  GT.  O)  CALL  SHIFTL<  X,  2000,  NP) 

CALL  FFOUR<  X,  1024,  C,  -1.  0) 

X MAX-CABS  <E<1)) 
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0055  DO  12  >1,512 

0056  X(l.  I>CABS<E<I>> 

0057  _ 12 _ IFiXHAX.  LX  X<1,  1))  XMAX»X<1,  I)  . 

0058  XS-1690.  O/XMAX 

0059  N*4090.  0-X  <1,1  )*XS 

0060  CALL  IXYPT<N,0,0) 

0061  00  13  >1,512 

0062  IXX«8#I-4 

0063 _ JJCMGS&  JJfcJLOjJl  #XS _ 

0064  13  CALL  IXYPT<  XY,  1XX,  1 ) 

0065  CALL  1 X  YPT  <  0,  4095,  0 ) 

0066  CALL  JXYPT <0,  0,  0) 

0067  CALL  IXYPT <  4095,  0,  0) 

0068  CALL  IXYPT <4095,  4095,  0) 

0069  _  GO  TO  1  _  _  _ _ 

0070  10  CALL  XY END 

0071  CALL  EXIT 

0072  END 


ROUTINES  CALLED: 

X1INIT,  SIZE  < _  CHSIZE,  IXYPT  ,  XYCHAR,  SMIFTL,  FFOUR 
CABS  ,  XYEND  .  EXIT 


SWITCHES  -  /ON 

BLOCK  LENGTH 

MAIN.  . .13422  .<  064334  >  * 

♦♦COMPILER  -  CORE*# 

PHASE  USED  FREE 

DECLARATIVES  00366  05535 
EXECUTABLES  00937  04964 
ASSEMBLY  01478  07340 


Attachment  #2  Application  Note  on  OCT 

DCT  -  Disci  be  Chebyshev  Transform 

A.  Implement  at  1-:-'  • 

The  Discrete  Chebyshev  Transform  of  a  data  set  y(m)  is  defined 

by  [1] 


TOO  -  £  R.  |  .-'“W 

N-l 

z 

} 

j 

meo 

k  -  0,1,...  ,N— 1 

end  A  ■  W  2 

k  ■  0 

1  2 

k  ■  1,2,..., N-l 

The  DCT  program  evaluates  this  expression  utilizing  an  in-place, 
UH-point  Past  Fourier  Transform;  the  various  parameters  employed  are 
as  follovs: 

N  *  number  of  data  pts  y(m) 

*  2P  where  p  la  Integral 

1(2  *  2N 

NU  -  UN 

HM  ■  number  of  sub-matrices  employed  in  the  FFT 

•  p  ♦  2 

Sign  ■  -1.0  for  forward  transform 

■  +1.0  for  inverse  transform 

Because  the  computations  are  done  in-place,  Y  is  both  input  and  output 
arrays;  X  is  a  complex  buffer  array. 

The  program  generates  2N  output  coefficients,  of  which  the  first  N 
contain  the  correctly-ordered,  desired  output.  Although  this  program 
utilises  an  FFT,  unlike  the  FFT  only  the  2N  real  output  points  from  the 
forward  transform  are  required  to  generate  the  Inverse. 

B.  Program  Notes 

The  FFT  employed  is  Robinson's  decimation-in- frequency  algorithmic  J. 

The  *phase-shift '  operates  on  the  complex  X  array  only  when  its  components 
represent  functions  of  k;  i.e.,  for  the  forward  transform  the  FFT  is  done 
first,  then  the  array  is  phase -shifted;  for  the  inverse  transform  the  reverse 
procedure  is  followed. 

Each  time  the  FFT  is  performed,  the  bit-reversed  outputs  ave  unscrambled 
within  the  subroutine  to  produce  a  sequentially-ordered  array;  the  output 
consists  of  the  real  component  of  the  complex  array. 
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DIMENSION  Z ( 1024) 

DIMENSION  Y < 1 024 ) <  MM<12) 

REAL  X  <2,2046 ),  B  ( 1 00 ) _ 

BYTE  IA(  70),  IX<9),  IQ<10),  TEST 
EQUIVALENCE  <E,  X) 


! 

j 


6 

7 

1 

20 

U 


_ DATA  TEST. I X , IQ/  X  ,  IE ' « 
♦  'E',  'A', 

CALL  XYINIT <B, 100) 

CALL  CHSIZE<0.  4.  90.  0 ) 
FORMAT  < 13) 

IC-0 


X'  r  'P'  f  O'l/S'  t/X  j  'O' i l  ' N I ' , 

'R'.  'T'»  'H'  »  O',  'U'»  'A',  'K' ,  'E'/ 


READ <6, 6) IF 

DO  20  1*1 ,1024 
Y  ( I )  *0.  0 


IC-IC+ 1  _ 

READ  <  I,  END-10)  I A 


READ<1)NPT,  <Y<I),  I-1,NPT) 
IF ( IF.  GE.  IC)  GO  TO  S 
REWIND  1 
GO  TO  7 


8 _ IF  (IF.  N£.  IC)  GO  TO  11 

CALL  IXYPT  <230,  0,  O) . 

CALL  XYCHARUA,  20) 

CALL  IXYPT <486,0,  0) 

IF<  IA<  70).  EG.  TEST)  GO  TO  2 
CALL  XYCHARUQ,  10) 

GO  TO  3 

2  CALL  X YCHAR  < I X, 9 ) 

3  CALL  IXYPT <300,0,  0) 

CALL  IXYPT <2190,  0, 1 ) 

DO  13  1*1, 1024 

15  Z< I )*Y< I ) 


— . „.£ALL  OCT <312,  1024,  2048,  11,  MM,  X,  Y,  -1.  O) 

XMAX=/<1) 

XMIN-XMAX 
DO  4  1*2,512 

IF  <  XMIN.  GT.  Y<  I ) )  XMIN-YU) 

4  IF<XMAX.  LT.  Y<  I ) )  XMAX-Y<I) 

.  XjS-1690.  0/  <  XMAX-XMIN) 

N-2190.  0+XMlN#XS 
CALL  I XYPT <N«  0, 0) 

CALL  IXYPT<N, 4000, 1) 

IXX-0 

IY-2190.  0+<XMIN-Y<l) >#XS 
CALL  IXYPTUY,  IXX,  0) 

DO  5  1-2,  NPT . ~ . "  . . . . . . 

IXX-4*I-4 

IY-2190.  0-MXMIN-Y<I))#XS 

5  CALL  IXYPT <  I Y,  IXX,  1) 

CALL  I  XYPT  <2390,  0,0) 

CALL  IXYPT < 4090,  0,  1 ) 

DO  9  I- 1,1024 
X<1,  I)«Z<I) 

X<2,  I)»0.  0 
DO  14  1-1025,2048 
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X  <  1*  I  )*0.  0 
14  X<2,  I)*0,  0 

_ *_162iLi _ 

XMAX»X ( 1,  1 ) 

XMIN-XM AX 
DO  12  1*1,  1024 

IF(XMIN.  GT.  X<  1#  2  )  )  XMIN»X(1,I) 

12  IF(XMAX.  LT.  X<  1,  I ) )  XMAX*X<1,  I) 

_ X  S«  169Q,  O/JJK  rjftX -A  niNJ _ 

N*4090.  0*XHIN#XS 
CALL  IXYPT<N,  0,  0) 

CALL  IXY?T(N, 4095/  1) 

N*4090.  0+<  XMIN-X < 1, 1))#XS 
CALL  IXYPT (N,  0,  0) 

_ DO  13  I«2,  1024 _ 

IXX«4#I-4 

IY-4090.  0+<  XMIN-X  (1,1)  )#XS 

13  CALL  IXYPT <  I Y, IXX, 1 ) 

CALL  IXYPT <0,  4095/  0) 

CALL  IXYPT <0,  0,0) 

_  CALL  IXYPT <4099,  0,0) 

CALL  IXYPT ( 4095,  4095,  0 ) 

GO  TO  1 

10  CALL  XYEND 

CALL  EXIT 
END 

ROUTINES  CALLED: 

XYINIT,  CHSIZE,  IXYPT  ,  XYCHAR,  DCT  ,  FWT  ,  XYEND 
EXIT 


SWITCHES  ■  /ON,  /SU 


BLOCK .  " ”  LENGTH 

MAIN.  13316  (064010)# 

♦♦COMPILER  -  CORE## 

PHASE  USED  FREE 

DECLARATIVES  00366  05529 
EXECUTABLES  00937  04958 
ASSEMBLY  01472  07340 


Attachment  #b 


Plot  of  DCT  and  WHT  for 
ACDA  Seismic  Signature  Data  Base 
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